home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / fortran / extenpr.zip / CALPI.FOR next >
Text File  |  1988-06-17  |  3KB  |  108 lines
  1.     PROGRAM CALPI
  2. C    This program calculates pi using a method described in 
  3. C    Scientific American February 1988, page 115.
  4. C
  5. C    By  Mark P. Esplin
  6. C        16 Shawsheen Rd.
  7. C        Billerica, Ma 01821    
  8. C
  9.     INTEGER*2 Y(3000),ALP(3000)
  10.     INTEGER*2 T1(3000),T2(3000),T3(3000),T4(3000),T5(3000)
  11.     INTEGER*2 T6(3000)
  12.     INTEGER*2 NT,N
  13.     INTEGER*2 ITER,NITER,ITEMP,IFIL
  14.     CHARACTER ANS*1,FILNAM*20
  15.     DATA NT/3000/
  16. 10    WRITE(*,*) 'Enter number of desired digits'
  17.     READ(*,*) NDIG
  18.     N=NDIG/4+2
  19.     IF(N.GT.NT) THEN
  20.       WRITE(*,*) 'Not enough memory available'
  21.       GOTO 10
  22.     ENDIF
  23. C    For now number of iterations is 
  24.     WRITE(*,*) 'Enter number of iterations'
  25.     READ(*,*) NITER
  26.     IFIL=0
  27.     WRITE(*,*) 'Log results to disk?'
  28.     READ(*,'(A)') ANS
  29.     IF(ANS.EQ.'y'.OR.ANS.EQ.'Y') THEN
  30.       IFIL=1
  31.       WRITE(*,*) 'File name?'
  32.       READ(*,'(A)') FILNAM
  33.       OPEN(1,FILE=FILNAM,STATUS='NEW')
  34.     ENDIF
  35. C    Set beginning values for Y and ALP.
  36. C    Y=SQRT(2)-1
  37.     ITEMP=2
  38.     CALL CVI2EX(ITEMP,T1,N)
  39. C    The square root use 4 temporaries, so it is actually using T3 - T6.
  40.     CALL EXSQRT(T1,T2,N,T3)
  41.     ITEMP=-1
  42.     CALL CVI2EX(ITEMP,T1,N)
  43.     CALL ADD(T1,T2,Y,N)
  44. C    ALP=6-4*SQRT(2)
  45.     ITEMP=4
  46.     CALL CVI2EX(ITEMP,T1,N)
  47.     CALL MUL(T1,T2,T3,N)
  48.     CALL NEG(T3,N)
  49.     ITEMP=6
  50.     CALL CVI2EX(ITEMP,T1,N)
  51.     CALL ADD(T1,T3,ALP,N)
  52. C    Main loop
  53.     DO 20 ITER=1,NITER
  54. C         T1=SQRT(SQRT(1-Y**4))
  55.       CALL MUL(Y,Y,T1,N)
  56.       CALL MUL(T1,T1,T2,N)
  57.       CALL NEG(T2,N)
  58.       ITEMP=1
  59.       CALL CVI2EX(ITEMP,T3,N)
  60.       CALL ADD(T2,T3,T1,N)
  61.       CALL EXSQRT(T1,T2,N,T3)
  62.       CALL EXSQRT(T2,T1,N,T3)    
  63. C         Y=(1-T1)/(1+T1)
  64.       CALL COPY(T1,T2,N)
  65.       CALL NEG(T2,N)
  66.       ITEMP=1
  67.       CALL CVI2EX(ITEMP,T3,N)
  68.       CALL ADD(T2,T3,T4,N)
  69.       CALL ADD(T1,T3,T2,N)
  70.       CALL COPY(T4,T1,N)
  71. C            INV using T4-T6
  72.       CALL INV(T2,T3,N,T4)
  73.       CALL MUL(T1,T3,Y,N)
  74. C         T1=1+Y
  75.       ITEMP=1
  76.       CALL CVI2EX(ITEMP,T2,N)
  77.       CALL ADD(T2,Y,T1,N)
  78. C        T2=ALP*(1+Y)**4
  79.       CALL MUL(T1,T1,T2,N)
  80.       CALL MUL(T2,T2,T3,N)
  81.       CALL MUL(ALP,T3,T2,N)
  82. C        T1=-(2**(2*(ITER-1)+3))*Y*(1+Y+Y**2)
  83.       CALL MUL(Y,Y,T3,N)
  84.       CALL ADD(T1,T3,T4,N)
  85.       CALL MUL(T4,Y,T5,N)
  86.       ITEMP=2**(2*(ITER-1)+3)
  87.       CALL CVI2EX(ITEMP,T6,N)
  88.       CALL MUL(T6,T5,T1,N)
  89.       CALL NEG(T1,N)
  90. C         ALP=T2+T1
  91.       CALL ADD(T2,T1,ALP,N)
  92. C         PI=1/ALP
  93.       CALL INV(ALP,T1,N,T2)
  94.       WRITE(*,*)
  95.       WRITE(*,900) ITER
  96. 900      FORMAT(' Iteration number ',I3)
  97.       IF(IFIL.EQ.1) THEN
  98.         WRITE(1,*)
  99.         WRITE(1,900) ITER
  100.         ITEMP=1
  101.         CALL WRTEXT(ITEMP,T1,N)
  102.       ELSE
  103.         CALL PRTEXT(T1,N)
  104.       ENDIF
  105. 20    CONTINUE
  106.     STOP
  107.     END
  108.